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ABSTRACT 

We present an analytical class of equilibrium solutions for the structure of relativistic 
sheared and rotating magnetized jets that contain no boundary current sheets. We 
demonstrate the overall dynamical stability of these solutions and, most importantly, 
a better numerical resistive stability than the commonly employed force-free struc- 
tures which inevitably require the presence of dissipative surface currents. The jet is 
volumetrically confined by the external pressure, with no pressure gradient on the 
surface. We calculate the expected observed properties of such jets. Given the sim- 
plicity of these solution we suggest them as useful initial conditions for relativistic jet 
simulations. 
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1 INTRODUCTION 

Jets are present in a great variety of astrophysical structures 
ranging from small scale protostellar jets ( [Reipurth et al.| 
|1998[ ) to large-scale jets feeding the enormous lobes in radio 
galaxies ( Rees||1978 1. In addition, there is a broad range in 



jet power and lifespan, for instance 7-ray bursts release most 



of their energy through a jet within a few seconds ( Sari et al 



19991, whereas quasar jets are long lived entities which last 



for millions of years. 

Measurements of synchrotron radiation and rotation 
measure suggest that the presence of magnetic fields is a 



ubiquitous element of jets (Gabuzda et al. 20041 and it is 



a standard ingredient of jet models i.e. Komissarov (19991; 



Leismann et al. ( 2005 1; Tchekhovskoy et al. ( 2008 K The sim- 



plest approach to a magnetic jet model is that of a force- free 
field. In this case the magnetic field dominates the jet. So, 
assuming ideal MHD, the magnetic field comes to an equi- 
librium within a few Alfven crossing times, and relaxes to a 
force-free state. However, any force-free magnetic field that 
is bound in space must have a current sheet on its boundary. 
A current sheet is an infinitesimal surface current which sep- 
arates the non-zero component of the magnetic field which 
is parallel to this boundary from the external medium. Al- 
though this is a viable state for ideal MHD and it is possi- 
ble to verify the stability of the field ( Woltjer||1958 1, if one 



takes into account dissipative effects, surface currents are 
critical (Taylor 19861. Physically, surface currents shall be 



either sources of instability or they shall dissipate making 
the transition from the jet to the external medium smoother. 
In addition their detailed study requires a micro-physical ap- 
proach upon which there is not a generally agreed picture, 
as in principle reconnection is described through different 
models ( |Sweet|1958| |Parker|1957} |Petschek|1964| [Uzdensky 
20TT|). 



A serious drawback of surface currents is their numer- 
ical treatment in simulations. A surface current formally is 
the derivative of a step function of the magnetic field, which 
is none other than a Dirac-5 function. Such discontinuities 
make this study a laborious task. In general, steep but fi- 
nite derivatives of the magnetic field are interpreted as sur- 
face currents. Observations of distant structures cannot pro- 
vide sufficient information about their presence. Solar sys- 
tem observations show rapid transitions in magnetic fields 
which are associated with large current densities, especially 
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in structures associated to coronal mass ejections (Burlaga 
^et al. 1981). Nevertheless, these are explosive entities, and 
if the dissipation timescale is longer than their dynamical 
evolution they can retain surface currents. 

The elimination of surface currents is feasible through 
the inclusion of plasma pressure in the dynamics of the prob- 
lem. In particular, when some gas pressure is included in 
the system the basic force-free equation is substituted by 
the Grad-Shafranov equation ( Shafranov||1966 1. This equa- 
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tion takes into account both the force of the magnetic field 
on the current and that arising from the pressure gradient; 
the equiUbrium state is given by the balance of these forces. 
Unlike the force-free description of the problem, this exten- 
sion allows a smooth transition from the magnetized area to 
the external medium which contains no magnetic field. This 
is consistent with the result of resistive decay, which turns 
magnetic field energy into heat, and thus, to an increase in 
pressure. 

It is also possible to study the relativistic generalisa- 
tion of force-free systems. Such studies have applications 
in the context of 7-ray bursts, AGN jets, and microquasar 
jets. In relativistic structures the force of the electric fields 
on the charges are significant and should be taken into ac- 
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progress has been made in this direction, in special and gen- 



eral relativistic formalism ( Komissarov 2002 , Gammie et al 
20031 IFendt fc Ouyed|[2004[ |McKi'nney||2006a|b| |McKinney 



fc NarayanpOOf l. Another field of application of force- free 



relativistic magnetic fields is that of pulsar magnctospheres 
(Contopoulos et al. 1999 Goodwin et al. [2004; S pitkovskyj 



20061. Finally, it is also possible to include pressure in such 
systems and take into account those dynamical effects which 
shall give the relativistic analogue of the Grad-Shafranov 
equation ( Gourgouliatos fc Vlahakis||2010 1. 

In the context of relativistic jet models for launching, 
acceleration, and coUimation, for years progress was made by 
relying on the assumption of steady state self-similar flows, 
e.g. Li et al.| ( 1992 I; Contopoulos (19951 or force-free fields 



(Fendt 19971). However, recent relativistic MHD jet forma- 



tion simulations have been capable of following the jet from a 
Keplerian disk to several thousand Schwarzschild radii above 
the disk (iPorth & Fendt||2010| iPorth et al.||2011|). Simula- 



tions in GRMHD resolve the accretion process towards the 
black hole, outflows, or the efficiency of a Blandford-Znajek 



tower, e.g. De Villiers et al. (20051; McKinney & Narayan 



( |2007| ). 

Unlike the above work which includes jet production; 
we focus on the asymptotic region of the jet and work in 



the context of magnetic towers ( Lynden-Bell 2003 Uzdensky 
fc MacFadyen|[2006 1 that do not have surface currents and 



where the field is confined inside a cylinder. Thus, our solu- 



tion differs from the approach of Appl & Gamenzind ( 1993 1 
who first presented a solution to the force-free asymptotic 
Grad-Shafranov equation that extends across the light cylin- 
der, but does not include return currents. As opposed to the 
magnetic tower case, in our work the field is not force-free 
but it coexists with a plasma whose pressure has some dy- 
namical contribution. This allows us to make the transition 
smoother from the area dominated by the magnetic field to 
the area dominated by the pressure. 'Gourgouliat os et al.| 
(20101 have found a similar class of solutions for fields of 



topology similar to a spheromak both for static and expand- 
ing structures ([Lyutikov fc Gourgouliatos||2011| IGourgorT] 
liatos fc Lyutikov|2011 1 . The solutions we have found satisfy 
the Grad-Shafranov equation with the additional constraint 
that there are not surface currents on the boundary. We ap- 
ply this idea to jets by solving the Grad-Shafranov equation 
in cylindrical geometry. We remark that the non-linearity 
of the Grad-Shafranov and the free parameters allow many 
possible equilibrium solutions, even if we assume some sym- 



metry. A possible approach is to assume that the detailed 
structure of the field is dictated by the behaviour of the 
source of the jet. This may be true in the region very close 
to the origin, but we expect that further out the jet will 
relax to a state where quantities such as the flux and the 
total helicity will be conserved without memory of the fine 
details of the jet launching region ( |Spruit|[20T0| ). The prob- 
lem at this stage is underdetermined, as we can construct a 
great number of solutions, since in principle we can choose 
a form for two of the three basic physical quantities appear- 
ing in the problem, i.e. poloidal and toroidal field, and solve 
for the third one. However any solution of this kind is not 
necessarily a physical configuration. A forward method of 
approaching the problem is the following. Consider a physi- 
cal system of a spinning disc on which the field is anchored 
and a pressure environment. The relation between and 
shall be dictated by the way the disc spins. Because of 
dissipation, fields without current sheets are more likely to 
occur provided there is sufficient time for the field to relax in 
a dissipative-structure. This demand leads to a system where 
both components of the field go to zero on the boundary. 
Then we can determine the pressure inside the jet through 
the Grad-Shafranov equation, but it also has to be such that 
the total energy carried by the system is a minimum for the 
given boundary conditions, so that the equilibrium is stable. 
Therefore the physical problem is summarized in the follow- 
ing steps. Some poloidal flux is connected to a disc, which 
spins and generates some toroidal flux. Inside the jet there 
is some pressure, which is lower in regions with a stronger 
magnetic fleld, and increases as one moves to regions of 
weaker magnetic field and dominates totally outside of the 
jet. The family of solutions of the Grad-Shafranov equation 
gives combinations of functions which correspond to equilib- 
ria, but not necessarily to the most economic ones in terms 
of energy. So we demand that the energy is also minimum 
for given boundary conditions. Although the above process 
is natural, it is far from being analytically soluble. Using 
the experience of force-free fields where stable solutions are 
given for fields where the same amount of current is carried 
by each field line, which are none others than the constant a 
solutions ( Tayler|1973 1, we shall relate linearly the pressure 
to the poloidal or the toroidal flux. Then, by performing 
some preliminary simulations we shall verify that the fields 
are not destroyed by instabilities. 

We consider these solutions useful for two reasons. First 
they may be realistic physical states after dynamical and dis- 
sipative relaxation has taken place, the dynamical relaxation 
leads to force equilibrium, while the dissipative relaxation 
leads to the elimination of surface currents and a possible 
scenario is that of mixing of the gas with the magnetic field. 
When the timescales of the jet are shorter than the ones of 
the external medium the system has time to relax into a 
stable equilibrium. This is true for a jet with lower density 
compared to the external medium, as the Alfven velocity, 
determining the timescale of the jet is faster than the sound 
speed of the denser external medium. Second we suggest 
that these magnetic structures can be used as trial solutions 
or limiting states in jet simulations. They are simple enough 
so that they can be used without significant modification in 
the existing simulations. In addition the simulations we have 
performed demonstrate that indeed, they are viable models. 

We shall present two types of systems. The first one 
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corresponds to the static problem where there is no mo- 
tion at all. Such configurations have only magnetic field and 
gas pressure which are in equilibrium. The second one is the 
treatment of relativistic outflows. The plasma moves parallel 
to the axis, but the time derivatives of the physical quanti- 
ties are zero, giving stationary solutions. This configuration 
contains both electric and magnetic fields and plasma pres- 
sure. The electric field is induced by motion parallel to the 
axis of the cylinder. Therefore the forces that come to equi- 
librium are the force of the magnetic field on the electric 
current, the force of the electric field on the charges and the 
gradient of the gas pressure. 



2 STATIC JETS 

In this section we set the mathematical framework of the 
non-relativistic static problem. The basic equation we solve 
is the Grad-Shafranov equation -j x B = Vp, the current 
is given by V x i? = — j. 



(V X B) X B = 47rVp, 



(1) 



where B is the magnetic field and p is the plasma pres- 
sure. The Grad-Shafranov differential equation is a non- 
linear partial differential equation and its solution is a rather 
complicated task. In general some assumption of symmetry 
is made, and the physical quantities appearing can be ex- 
pressed in terms of two coordinates, but still the field can 
have components in all three dimensions. The physical quan- 
tities appearing are the poloidal field, the toroidal field and 
the pressure. In this paper we examine the most fundamen- 
tal case, where the field is cylindrically symmetric, namely it 
is symmetric under translations and rotations with respect 
to an axis, therefore the physical quantities shall depend 
only in the radial coordinate. We shall express the field in 
two different ways. The field written in terms of the poloidal 
fiux Pp and poloidal current Ip is 



B = \/Pp{R) xV<t> + Ip{R)V(j} , 



(2) 



and we can also write it in terms of Pt which is related 
to the toroidal fiux: J B^dR = Pt(-R)-l-const., by an addi- 
tive constant, note that as the system is symmetric under 
translations in z there is no need to integrate along this 
coordinate. 



VPt(i?) X Vz + It{R)Vz. 



(3) 



In the absence of pressure the field comes to a force-free 
equlibrium which is given by V x B = aB. The most sta- 
ble fields are given for a spatially constant a ( Tayler|1973 1. 
When we write the fields in terms of the fiuxes we find 
that the constant a solutions correspond to fields for which 
Ip,t = ap,tPp,t- The final solution, independent of the choice 
of the representation is that of the [Lundqulst] ( |1951[ ) field: 

Bff = coa{Ji{aR)4) + Jo{aR)z) , (4) 

where co is a normalisation constant and Jo and Ji are the 
Bessel functions of zeroth and first order. Fig. ([T]). 

In the Grad-Shafranov equation there is an extra force 
arising from the pressure gradient. Analytical solutions are 
possible under the assumption of a pressure profile pro- 
portional to the fiux. Unlike the force-free case, the solu- 
tions we find are different if we take the pressure to be 




Figure 1. The Lundquist force-free magnetic field. The solid line 
is the component of the field and the dashed is the . There 
is not a natural end of the confined field, one can choose to end the 
jet where the B^ field becomes zero for the first time _R = 1, in this 
case the component has already reversed its polarity, or earlier 
where at R = 0.63 where the Bz becomes zero for the first time 
(vertical line). The fields never become zero simultaneously thus 
there is no bound magnetic field configuration without surface 
currents. 



associated with the poloidal or toroidal flux. So we write 
Pp.t = 4^-Fp,t-Pp,t+Po p,t, where Fp^t are constants, andpo p,t 
is an additive constant, which does not have any dynami- 
cal role in the problem, but it has a lower limit to ensure 
that the plasma pressure is positive everywhere. Substitut- 
ing those in equation ([TJ and solving for the two cases we 
find that 



Pp = CpRJi{apR) 



Ft 

Pt = ctJo{atR) 2 ' 



(5) 



(6) 



where Cp,t are normalisation constants. We are looking for 
solutions without surface currents confined to a cylinder 
of unit radius Rjet = 1, thus in both solutions we want 
Pp,t{l) = and Pp_t(l) = 0. Applying these conditions and 
choosing a normalisation so that the maximum value of the 



field is unity we find Cp — 0.172, ap 
and ct = 0.186, at = 3.83, Ft = -1.10. 
give two physically distinct fields 



5.14, R 



-1.54, 



{cpapJo{apR) 



2F 



These two solutions 



(7) 



Bt = CtOitJi{atR)<p ~\- {ctatJo{atR) )i , 

at 



(8) 



which are plotted in Figs. ([2j|3|). The value for plasma /? — 
Snp/B'^ that we take for the lowest value of po p,t is plotted 




Figure 2. The magnetic field arising from Pp, the solid line is 
the Bij, component of the field and the dashed is the . Note 
that Bz changes polarity before reaches zero, resembling the 
behaviour of the force-free Lundquist field. 




R 



Figure 3. The magnetic field arising from Pt, the solid line is the 
B^ component of the field and the dashed is the Bz- Note that 
these fields correspond to physically distinguishable structures 
from the case presented in Fig. |2| and the force-free Lundquist 
field, Fig. g. 

in Fig. Q. For both configurations, in the greatest part of 
the jet /3 < 1, including regions where /? << 1 where the jet 
is practically force-free. Even if we do not choose to express 
the field in terms of some fiux we can write B — B^(f> + BzZ 
and the Grad-Shafranov equation is written in the form 

^(B^z+BI + SttpY + Bl = 0. (9) 




Figure 4. Logarithmic plot for the value of plasma /3 for the 
lowest acceptable value of po a-s a function of R in units of 
Rjct- The solid red curve corresponds to the poloidal solution 
and the blue dashed curve to the toroidal solution. The value of 
/3 becomes zero at 0.52 Rjot for the poloidal-type solution; and 
on the axis for the toroidal-type solution as the plasma pressure 
is zero there. In contrast, it becomes infinite as i? — > Rjct where 
the magnetic field is zero. 

Despite the fact that this equation does not add new infor- 
mation in the problem it shall give useful insight when we 
approach the outfiow problem. 



3 SHEARED JETS 
3.1 General solution 

In this section we study outfiows parallel to the axis. This 
procedure can work for any velocity, from slow ones to rela- 
tivistic. To model the problem, we keep the cylindrical con- 
figuration of the static problem but we impose a velocity 
along the z-axis v = VzZ, whose profile is cylindrically sym- 
metric, independent of time, but varies with radius. This 
velocity induces a radial electric field, whose divergence cre- 
ates an electric charge density and electric forces appear. 
As the magnetic field has no R component the field lines 
are not stretched and retain their shape. In addition there 
is no time dependence thus we do not expect displacement 
currents. The magnetic field is 

B = B^{R)4> + BziR)z, (10) 
and the electric field 

E = Er{R)R, (11) 
which is associated to the velocity by 

Er = VzB^ . (12) 

The time derivative of the magnetic field is zero and the 
electric field is curl-free, thus the Faraday induction equation 
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\7 X E- 



1 SB 

c dt 



holds. The electric charge density is 
E, (13) 
the electric current density is 



•0 V7 



(14) 



and we can verify that Gauss's law for the magnetic field 
V ■ B — holds. The momentum equation is 



rpo(^ + ^^- v)(c^1;)-Vp + 



/'i^- 



x B 



= 0, 



(15) 



where F = (1 — \)~^''^ is the Lorentz factor, po is the rest 
mass density and ^ is the relativistic specific e ntha lpy over 
c^. The first term of the momentum equation (15 1 appear- 



ing in brackets is zero as there is no time variation and the 
velocity is along z so from the V operator we shall keep only 
the derivative with respect to z, but since our system does 
not have a z dependence its contribution is zero. We remark 
that if we chose another velocity profile which contained a 
<j) component, inertial forces should appear. This case is ex- 
amined in greater detail in the simulation part of the paper. 
In the analytic part we focus on the axial motion, where the 
equilibrium is reached through the balance of field forces and 
pressure. For the present we do not need to take into account 
the details of the mass density and the relativistic specific 
enthalpy of the system provided that they do not vary with 
time or z. Thus we are left with the last two terms involving 
the pressure and the electromagnetic fields. Now we proceed 



to the solution of equation ( 15 I which after substituting the 
fields can be written as 



El + Svrp)' + Bl-El 



0. 



(16) 



Note the similarity with equation ([9|, with the extra term 
—E\ which is a result of the relativistic invariance of — 
E^ . For that reason, we shall introduce the new variable 



(17) 



which can be substituted to equation ( |16[ ) and it shall give 
R, 



-{Bt + + 8TTpy + = 0. 



(18) 



Based on the results of the previous section, we can gener- 
alise the solutions we have found in the previous section. In 
the static regime we had two representations corresponding 
to the flux functions Pp and Pt . Here instead we shall use the 
generalised functions Gp and Gt- The first representation is 



H 



R 
1 dGp 



R dR 



— FG 



'Pp,0 ■ 

The second representation is 
dGt 



(19) 



H 



47r 



dR ' 

B, = atGt , 

FtGt+pt,o- 



(20) 



The solution to the differential equation is exactly the same 
as the one in the static problem and is given by equations 



^ and ([gJi, while the values of Fp^t and ctp.t are the ones 
discussed there, for the same boundary conditions, namely 
the boundary of the jet is at i?jct = 1 and does not have sur- 
face currents. Then we can evaluate the electric field and the 
azimuthal component of the magnetic field, from equations 



(121 and (17 1 



Bj, 



El = 



(21) 



(22) 



Now the issue is the velocity profile we choose which gives 
the partition of H in the electric and magnetic field. 



3.2 Solutions for a given velocity profile 

Having solved the problem for the generic velocity profile 
we can write the expressions for the fields. The field and the 
pressure for the poloidal-flux approach are: 



Er 



Ba, 



CpapJi{apR) ' 
'"p 



(l-«2)l/2 



CpapJi{otpR) 



FpR\ 



Bz = CpUpJo^apR) 



2K, 



1 / F \ 
P= -^Fp[cpRJi{apR) +Pp,o, 

and for the toroidal-flux approach are: 



-Efl = -jz \rTi^CtatJi{atR) : 

(1 - viy'^ 



B. 



1 



(i-«i)i/ 

Bz = at (^ctJo{atR) 



Yj^ctatJi{atR) , 



^,2 I ' 



1 / F \ 



(23) 
(24) 
(25) 
(26) 

(27) 
(28) 
(29) 
(30) 



For instance, we choose a linear velocity profile which is 
maximum on the axis and zero at _R = 1, = Uz,o(l — ^)- 
For this particular case the fields have the forms plotted in 
Figs. ([5| and (|6|. Note that the values of Qp,i and Fp^t are 
those found in the static problem, so that the jet is confined 
to a cylinder of unit radius with no surface currents. 

We remark that v is the velocity of the outflow, and is 
not the same as the velocity of the field lines. The velocity 
of the field lines in units of c is given by 



ExB (-V X B) X B V B^ 

^^ = -^- W ="-T^' 

thus for V along the z axis it is 



Vf = 



+ 



VzBi 



(31) 



(32) 



B2 ' 52 

There is a difference between v and Vf by the projection 
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B, E,p 





Figure 5. The magnetic and electric fields and the pressure aris- 
ing from Gp. The solid line is the -B^ component of the field, the 
dashed is the Bz, the dotted the electric field and the dotted- 
dashed line is the p, on which a suitable value for pp^Q has been 
chosen so that it is positive, note that the fields and the pressure 
are of different dimensions so there shall be no comparison in the 
plot. 



B, E, p 




Figure 6. The magnetic and electric fields and the pressure aris- 
ing from Gf. The solid line is the S0 component of the field, the 
dashed line is the Bz, the dotted line is the electric field and the 
dotted-dashed line is the p, on which a suitable value for pt,o has 
been chosen so that it is positive. Note that the fields and the 
pressure are of different dimensions so there shall be no compar- 
ison in the plot. 



Figure 7. The velocities of the flow and the field lines normalised 
to c for the field found in the poloidal-flux approach. Fig. The 
solid line is the toroidal component of Vp, the dotted line is the 
axial component of Vp and the dashed line is v which is chosen 
to be in the z direction. In this case we have chosen that the 
maximum velocity for v to be 0.5c. 



of V on B. Vp has both an axial component and a toroidal 
component. Figs (|8|. Note that the plasma can drift 
along the field lines, thus it is not necessary that it moves 
with the Vp velocity, but in principle it may move with any 
velocity which is the sum of Vp and a component parallel 
to the magnetic field. The addition of a velocity component 
parallel to the magnetic field has no effect in the dynamics 
of the problem in the absence of inertia or for an appropriate 
drift along the field lines so that the particles perform no az- 
imuthal motion ( Gourgouliatos fc Lynden-Bell|20"TT |. When 
inertia is taken into account the dynamics of the problem 



change as the first term in equation ( 15 1 gives a centrifugal 



force when there is an azimuthal component in the velocity. 
We consider this force in the simulations below. Note also 
that in the demand of a smooth transition from the jet to 
the external medium Vz needs to go to zero at the boundary 
of the cylinder, otherwise there will be a discontinuity in 
velocities in this layer. 



4 MHD SIMULATIONS OF JET 
PROPAGATION 

One of the main advantages of the magnetic field distribu- 
tion proposed in this paper is that it provides a transition 
through the jet boundary i?jct = 1 where no surface force is 
present. 

Simulations of jet propagation for helical magnetic fields 
often rely on a toroidal magnetic field distribution which in- 
troduces spurious relaxation processes along the jet bound- 
ary. The problem lies in the fact that one usually injects a 
magnetized jet with a helical field into an ambient gas ini- 
tially without any toroidal field. Examples of this approach 



are the seminal papers by Clarke et al. ( 1986 1 and Lind et al 
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Figure 8. The velocities of tlie flow and the field lines normalised 
to c, for the field found in the toroidal-fiux approach, Fig. (j6]l. The 
solid line is the toroidal component of Vp, the dotted line is the 
axial component of Vp and the dashed line is v which is chosen 
to be in the z direction. Note that we use the same profile for i; 
as in the case of Fig. ^ where its maximum value is 0.5c. 



(|1989|), but also many follow up simulations ([Koessl et al 
1990 ; Fr ank et al. (19981 |Stone fc Hardee|2000| [o'Neill et at 
2005)), and more recently the relativistic cases presented by 
Keppens et al.| ( |2008| ); [Mignone et al.| ( |2010| ). Other au- 
thors have applied a force-free field distribution over the 
whole computational domain ( Todo et al.|1992 19931, i.e. a 
toroidal field also in the ambient gas. 

For example, it is possible to choose a sinusoidal func- 
tion which leads to a vanishing toroidal magnetic field at the 
jet boundary, but this is not a force-free field configuration. 
Thus, Lorentz forces within the jet will distort the injected 
jet material. Furthermore, applying a vanishing toroidal field 
outside the jet (for R > Rjet), the derivative dB^/dR be- 
comes infinite and thus the Lorentz force at this radius. An- 
other option is to choose a linearly increasing toroidal field, 
which is force-free across the jet and maintains the injected 
jet structure, however, since it is set to zero for R > Rjct, it 
introduces again an infinite Lorentz force. 

A more elaborate initial and boundary condition has 



been applied by Keppens et al. (20081, who solved the full 



MHD equations to find a jet inlet in radial force-balance. 
However, even in this case the toroidal field distribution is 
cut off at the jet radius, potentially leading to relaxation 
forces along the jet boundary. 

In this section we test the behaviour of the magnetic 
field configuration derived earlier in this paper in simula- 
tions. We have therefore performed axisymmetric relativistic 
MHD simulations of jet propagation applying the PLUTO 
3.01 code ( [Mignone et al.|[2007) . We inject a jet with the 
velocity, pressure, and magnetic field profiles as derived in 
Section |3.1| into an ambient corona of vanishing magnetic 
field and uniform pressure. In addition to the properties 
derived, we need to define two other dynamical parame- 



ters which are the jet density and the density of the am- 
bient gas, Pjct,Poxt. We do this by applying a polytropic 
equation of state, thus obtaining Pjct(-R) = Pjct{Ry^'' / Kjct, 
and pcxt{R) = Pext(-R)^''"'/A"cxt. Note that depending on the 
choice of Kjet, Kex.t there will be an entropy jump between 
the jet and the external material. Essentially, we inject jets 
with internal Alfven Mach numbers of 2 — 10 and fast mag- 
netosonic Mach numbers of > 1, as expected from MHD 
disk jets. Clearly, the hydrodynamical parameters will play 
an essential role in the stability characteristics of these jets. 
However, a full stability analysis is beyond the scope of this 
paper as it would require also a 3D-treatment. Here we re- 
strict ourselves to a preliminary study of the morphological 
stability of the jet-ambient medium interface. For our sim- 
ulations we apply a numerical grid of physical size (5 x 20) 
with a resolution of 100 equidistant cells between r = and 
r = 1.5 and 200 scaled cells between r = 1.5 and r = 5. In 
2-direction we have 100 equidistant cells till 2 = 3 and 250 
scaled cells till z = 20. We also have performed comparison 
simulations with double resolution which did not show any 
significant differences. 

The parameters of our simulation runs are listed in 
Tab. [l] For comparison we provide as derived parameters 
the maximum poloidal Alfven Mach number 



Ma = 



^7. 



and the maximum poloidal fast "cold" magnetosonic Mach 
number 



Mf 



47rpjotuf 



with the relativistic velocity defined as u = Tv and the 
Lorentz factor P. 

To demonstrate the overall force-balance across the jet 
boundary, we plot the radial profile of radial Lorentz force 
component in Fig. [9] Inside the jet the Lorentz force is non- 
vanishing, however it is balanced by the jet kinematic forces. 
Towards the jet boundary the radial Lorentz force vanishes, 
and remains zero outside the jet. 

Figure [To] shows the dynamical parameters for our ref- 
erence simulation S06 at 500 dynamical time steps. The 
density and velocity distributions clearly demonstrate the 
smooth transition across the jet boundary, thanks to the new 
analytical MHD solution introduced above. Along the axis 
refiection shocks develop, which, however, do not strongly 
affect the jet/ambient gas interface. 

Figure [TT] shows the poloidal Alfven Mach numbers for 
simulations S04, S05, and S07. This figure displays a smaller 
subset of the computational domain. The fiows are all well 
super Alfvenic, while only a small amount of the entrained 
material along the fiow is sub Alfvenic. 

The appearance of refiection shocks is a well-known fea- 



ture in jet propagation simulations (see e.g. Hardee & Nor- 
( 1988) [T989|; [Hardee et al.[([l992|);[Bodo et al.[([2004[)) 



Following Hardee & Norman ( 1988 1989 1 we estimate the 



wavelength A of the reflections as 



A : 



Mj,in 



zR 



(33) 



m + 1/2 1 + ^ 
with the jet internal Mach number Afj.int, the density con- 
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o 



a; 




log(rho(r,z)) 



Figure 9. Radial profile of Lorentz force radial component for 
simulation S07 after 500 dynamical time steps, indicating a 
smooth transition from jet to ambient medium. Shown is the pro- 
file along the first active zones. 



trast r\ = pjet/pe^i, the wave mode number m, and half- 
thickness radius R. With the usual definition of the gener- 
alized relativistic Mach number M = Tjetv/TaVB with the 
Lorentz factors for the jet bulk speed Fjot = (1 — vfct)~^^^ 
and the sound speed Fs = (1 — I's)^^''^, and Vs is the sound 
speed normalized to the speed of light = a/c, and with 
the sound speed 

1/2 

(34) 



7P 



P + iP/il - ] 

I Hardee et al.|200T | For simulation run S07 we find 77 = 0.5, 
Pjet — 1.48, a sound speed a ~ 0.19, and with the Lorentz 
factors Fs ~ 1.02, Fjet — 1.25, resulting in a generalized 
Mach number M ~ 3.7, and, thus, a zeroth mode wave 
length A ~ 77?. If we take R = 0.5, the radius where we mea- 
sured the above-mentioned quantities, the predicted wave- 
length is A ~ 7R, which is similar to our simulations result 
(see Fig. 10). 

Figure [12] demonstrates how the jet penetrates the am- 
bient medium for simulation S04, showing the log-scale den- 
sity distribution of the whole computational domain at time 
200 and 500, respectively. This is an over dense jet with high 
(maximum) Lorentz factor F = 7.1 (i.e. iimax = 0.99). Again 
we see a smooth transition between jet and ambient gas with 
only little entrainment going on. In the case of this powerful 
jet, the wavelength of the reflection shocks is longer. 



5 MHD STABILITY 

The stability of MHD equilibria has received much atten- 
tion in the plasma confinement literature due to the violent 
instabilities to which Tokamak devices are prone (see e.g. 
|Kadomtsev|p66l |Bateman||T978| [FreidberglpST] ). A com- 
mon technique for evaluating stability involves introducing 
perturbations of the form oc eiq)i(m(f) + kz — uit) on top of 
the equilibrium configuration to examine how this affects 
the total MHD energy. Analyses of this type have produced 
two well-known analytical criteria that provide simple tests 
for stability, the Kruskal-Shafranov (KS) criterion and the 
Suydam criterion. 
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Figure 12. Axisymmetric jet propagation simulation S04. Large- 
scale density distribution after 200, 500 dynamical time scales, 
respectively. 
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Figure 10. Axisymmetric jet propagation simulation S06. Density p, axial velocity Vz, and relativistic axial velocity Uz = Tvz distribution 
after 500 dynamical time scales. Note the size of the figure which is a 2 X 10 subset of the whole 5 X 30 computational domain. 



The KS criterion is a test for the current-driven |m| — 
1 mode instability. This so called "kink" mode is one of 
the most dangerous instabilities for magnetically dominated 
astrophysical jets because, unlike other modes, it displaces 
the fluid center of mass by distorting the jet cylinder into a 



helix-like configuration, thereby converting magnetic energy 
into kinetic energy. The internal and external versions of 
this instability have already been the object of a number of 



studies in the astroph ysical literature (e.g.|Istomin fc Pariev 
[19961 [BegelmanI 19981 [Nakamura & Meier|2004| [GianniosX 
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Figure 11. Axisymmetric jet propagation simulations S04, SOS, S06, S07. Shown is the poloidal Alfven Mach number M\{r,z). Note 
the (real) rather long wavelength of the reflections, as only a subset of the whole computational domain is shown. 



Spruit||2006l [Narayan et al.||2009| ). The KS criterion states 
that the jet is kink stable if 



2nR B 
Li B, 



> 1, 



(35) 



where Lj is the length of the jet. [Narayan et aL| ( |2009[ ) have 
found that this criterion is similar in form for force-free rela- 
tivistic jets as long as one uses the jet frame magnetic fields 
in the above expression. Since realistic jets are conically 
expanding, ideal MHD suggests that B'^/B^ <^ 1, render- 



ing the kink mode particularly dangerous to jet stability. 
However, we note that a dynamically significant amount of 
plasma pressure makes the jet less susceptible to the kink 
instability. 



The Suydam criterion concerns the specific interplay be- 
tween the curvature of the magnetic field lines and the pres- 
sure gradient. For the stability of these local "interchange" 
modes, it is necessary (but not sufficient) for the following 
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Table 1. Relativistic MHD simulations of jet propagation, ap- 
plying model Fig. 5. Parameters common for all simulations are 
Fp = -1.54, Op = 5.14, Cp = 0.172, Ft = -1.10, at = 3.83, 
ct = 0.186. Parameters specific for each simulation are listed 
below, poxt = Pp,o or Poxt = Pt,o, respectively, tijet = fz,o- 
i^ext = P(Jxt- ^jct = Pjcl^ ^'^^ comparison we provide as derived 
parameters the internal Alfven Mach number M\, and internal 
fast magnetosonic Mach number MpMi given as a typical value 
along the jet (naturally these parameters vary substantially across 
the jet according to the velocity and magnetic field profiles, see 
Fig. 10). 



if M > /3 such modes are stabilized, though other global 
modes are excited with slow growth rates (for a discussion 
of pressure-driven modes in the context of jets, see |Lon-| 
[garetti_^2008| ). Thus, the apparent stability of the poloidal 
equihbrium in our jet propagation simulations may be due 



to high velocity shear, though Bondeson et al. ( 1987 1 's work 



would need to be extended to the relativistic regime for a 
more thorough understanding of how velocity shear affects 
pressure-driven instabilities in relativistic jets. 
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equality to be satisfied for the entire jet profile: 



RBt 
An 



(I) +8p'>0. 



(36) 



where q = ^ ■ Equation ( 36 1 shows that magnetic con- 
figurations with high enough values of q'^ can overcome 
the destabilizing effects of a negative pressure gradient (i.e. 
p' < 0) . Evaluating equation ( 36 1 for both equilibria reveals 



that the poloidal equilibrium is unstable between 7? = to 
0.4, and the toroidal equilibrium is stable for all R. 

As suggested by the stability of our solutions in the 
MHD propagation simulations (see the results from the 
KS and Suydam criterion should be approached with some 
skepticism. The kink mode, mostly discussed in the AGN 
jet literature in the context of magnetically dominated jets, 
is not as important in jets where the plasma pressure is dy- 
namically significant, as is the case with our solution where 
the plasma /3 (= Sivp/B^) is of order unity in the outer re- 
gion of the jet. In the case of magnetically dominated jets, 
simple analytical kink mode analyses which suggest they are 
kink unstable (e.g. |Begelman|1998[ ) do not take into account 
complicating factors such as causality in conically expand- 
ing jets, gradual velocity shear, and field line rotation. This 
may explain why some analytical stability analyses suggest 
jets are kink unstable while, in contrast, three-dimensional 
numerical simulations which do take into these complicating 
factors can produce jets which are kink stable ( |McKinney| 
|& Blandford 2009). In the case of pressure driven instabili- 



ties, the poloidal equilibrium is unstable per a naive appli- 
cation of the Suydam criterion. Unfortunately, this criterion 
does not take into account bulk velocity and therefore has 
limited applicability for jets. In the fusion confinement lit- 
erature, [Bondeson et al.| ( |1987[ ) have extended the Suydam 
analysis to include velocity and have shown that the stability 
of localized pressure-driven modes depends on the quantity 
M = v'^p^^ {B^q' /qy^ , a form of Alfvenic Mach number. If 
M < P then velocity shear destabilizes local resonant modes; 



6 SYNCHROTRON JET PREDICTIONS 

The polarization and anisotropy of the resulting synchrotron 
emission from magnetic jets containing relativistic plasmas, 
allows high resolution telescopes to probe the structure of 
jet magnetic fields. To do this, the transverse structure of 
the jet must be resolved, a feat that only very long base- 
line interferometric (VLBI) radio telescopes can perform on 
some active galactic nuclei jets. The emission patterns of a 
transverse (compared to the jet symmetry axis) cut on a jet 
threaded with a large-scale helical magnetic field has already 
been studied using analytical models (e.g. Laing|1980 Lyu- 



|tikov et al.|2005 Clausen-Brown et al.|20lT I and numerical 
simulations (e.g. Broderick & McKinneyj|2010j [Forth et al . 
|2011). Following Clausen-Brown et al 



72011) 



we create po- 
larization profiles by calculating the Stokes Q and I {U = 
since our chosen reference position angle is the jet direction 
projected onto the sky). We assume the entire jet cylinder 
is filled with emitting electrons with a power-law electron 
distribution function, dn — K^E~^dE, where the density 
of particles does not vary with R. Retaining the sign of Q 
in calculating 11 gives the synchrotron radiation's electric 
vector position angle (EVPA) as well: when Q is positive 
the EVPA is along the projected jet direction and when Q 
is negative the EVPA is perpendicular to the projected jet 
direction. Faraday rotation measure (RM) profiles are calcu- 
lated as well using a Faraday rotating electron density that 
is constant with radius. 

Figure [13] shows the calculated polarization and 
RM profiles using our solutions (labeled "Poloidal" 
and "Toroidal" ) and the Lundquist field with two different 
cutoffs chosen ("Force- free I" and "Force-free II") so they 
most resemble our solutions. The similarity between our so- 
lutions and their force-free counterparts will make distin- 
guishing them using VLBI profile difficult, especially con- 
sidering that most AGN jets are barely resolved if at all. 



7 DISCUSSION 

In this work we have considered the magnetic field structure 
of a magnetized cylindrically symmetric jet. We have found 
analytical solutions for the field and the plasma pressure 
that is confined by an external medium without a singular 
current sheet, but through a smooth transition. The configu- 
ration satisfies the Grad-Shafranov equation subject to given 
boundary conditions, namely both the magnetic fiux and its 
derivative have to be zero. Our assumption of cylindrical 
symmetry is valid in a jet's asymptotic region well down- 
stream of the coUimation region where the structure transi- 
tions from a broad outfiow into a conical one with a small 
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Figure 13. In the left column are plots of the fractional polar- 
ization of the jets as a function of projected radius for jet viewing 
angles of Sq;, = l/(2r). The sign of the polarization refers to 
whether the associated EVPA is parallel (Q/I > 0) or perpen- 
dicular {Q/I < 0) to the jet axis. The right column are plots 
of the Faraday RM that are normalized by the maximum change 
{ARM = RMmax—RMmin) in RM. Each of the four rows of plots 
refers to (a) the Lundquist solution where B^{R = 1) = 0, (b) the 
Bp equilibrium, (c) the Lundquist solution where Bz{R = 1) = 0, 
and (d) the Bt equilibrium. Notice that the magnetic field struc- 
ture for solutions (a) and (b) are similar in that they undergo 
a Bz reversal, and solutions (c) and (d) are similar in that no 
magnetic field reversals take place. It is for this reason that the 
plots on the first two rows resemble each other and the third and 
fourth row plots resemble each other. 



opening angle. Indeed, a jet powered by a magnetic mecha- 
nism ([Blandford fc Znajek|[T977! 'Blandfo rd fc Payne|[T982{ 
Lynden-Bell|2003[|Uzdensky fc MacFadyen|2006| ) shall start 
with a weakly collimated conical shape and, further down- 
stream, become c ollimated ([Vlahakis fc Konigl|2003||Komis- 



sarov et al.|2007| [Forth fc Fendt|2010| ). This is supported by 
observations of highly collimated jets observed in broad va- 
riety of astrophysical structures ranging from Herbig-Haro 



objects to AGN jets (Livio 1999 Reipurth fc Bally 20011 



M87 contains a well studied example of a collimated jet in 
which [Junor et al.| ( |1999[ ) have observed a broad structure 
very close to the origin with a half opening angle of 60° ; col- 
limation starts at around 30-100 Schwarzschild radii from 
the black hole and continues out to 1000s of Schwarzschild 
radii. The opening angle of the M87 jet at the parsec scale is 
around 6 — 7° indicating recollimation ( Kovalev et al.|2007 |. 
[Pushkarev et al.| ( |2009[ ) have found from a statistical study of 
relativistic jets that their intrinsic opening angles are around 
2 — 3°. Obviously, a global jet model shall account for these 
variations, but the use of a cylindrical geometry locally for 
regions where the opening angle is small, is a reasonable 
approximation. In our solution we have imposed a constant 
pressure environment which confines the jet to a cylinder. A 
more realistic model shall contain a varying pressure which 
decreases with distance from the origin and leads to jets of 
more complicated profiles, i.e. conical or paraboloidal. From 
an other point of view, close to the launching region the 
physical processes do not simplify to force equilibrium as we 
have assumed in our model, making our solution applicable 
away from the launching region. 



Within the limits of our approach our solutions might 
be interesting in the hght of the FR-I/FR-II dichotomy. 
FR-I sources have wider opening angles and entrainment 
with their environment and are therefore more unstable 
to Kelvin-Helmoltz modes (e.g. Bicknell|1995 I, while FR-II 
sources are more collimated and have a smoother interface 
with their environment in a manner more consistent with 
our solutions. 

Our solutions retain the choice of constant a, usually 
applied in force-free structures, because it is known to give 
the most stable amongst the force-free fields. The linear re- 
lation of the pressure on the fluxes allows analytical solu- 
tions. This gives two classes of fields: one where the axial 
field is reversing and another where the axial field always 
points in the same direction. These solutions can be gen- 
eralized to include relativistic outflows. We then use these 
solutions as starting points in simulations. Through these 
preliminary simulations we verify their stability and we also 
check the effects of centrifugal forces, as the particles do not 
move parallel to axis but follow the helical motion of the 
field lines, leading to a slightly different profile compared 
to the analytical solution. Finally the observational profiles 
predicted by this structure give reasonable results, but not 
clearly distinguishable from the ones predicted by force-free 
structures. 

We consider these solutions as a guide to investigate 
realistic configurations. They are simple enough and ana- 
lytical so that they can be used to parametrize physical 
systems or applied as trial cases in numerical simulations. 
Given an appropriate choice of initial magnetic field, mass 
load and velocity profile a jet can be simulated. Possible ap- 
plications are relativistic jets emanating from quasars and 
microquasars and in the non-relativistic context protostellar 
jets. 
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